Objetivo: aprender a manipular conjuntos de datos grillados y datos vectoriales mediante la aplicación de funciones en R.
Los datos vectoriales almacenan la ubicación, forma y atributos de características geográficas en una de las siguientes formas:
Existen diferentes paquetes que se pueden utilizar para cargar y manipular datos vectoriales (shapefiles) en R. En esta presentación, utilizaremos los paquetes sf y terra. El paquete tibble también puede ser instalado ya que es utilizado por sf.
Para instalar estos paquetes podemos utilizar la función install.packages.
Una vez instalados, podemos cargarlos con la función library.
Para leer datos vectoriales en R, podemos utilizar ambos paquetes. Para utilizar el paquete sf podemos aplicar la función read_sf.
Si la carpeta contiene más de un shapefile, hay que especificar la ruta del archivo:
Si escribimos el nombre del objeto en la consola:
countries
## Simple feature collection with 15 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS: WGS 84
## # A tibble: 15 × 10
## STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng
## <chr> <chr> <int> <chr> <int> <int> <dbl>
## 1 Member State NO 205 Rwanda 1000 3000 8.13
## 2 Member State NO 133 Kenya 1000 3000 48.5
## 3 Member State NO 74 South Su… 2011 3000 46.9
## 4 Member State NO 257 United R… 1000 3000 83.0
## 5 Member State NO 253 Uganda 1000 3000 23.2
## 6 Sovereignty uns… YES 61013 Ilemi tr… 1000 3000 2.84
## 7 Member State NO 79 Ethiopia 1000 3000 50.4
## 8 Member State NO 77 Eritrea 1000 3000 50.0
## 9 Member State NO 43 Burundi 1000 3000 9.12
## 10 Member State NO 68 Democrat… 1000 3000 99.0
## 11 Sovereignty uns… YES 40760 Hala'ib … 1000 3000 8.51
## 12 Sovereignty uns… YES 40762 Ma'tan a… 1000 3000 2.06
## 13 Member State NO 6 Sudan 2011 3000 81.9
## 14 Member State NO 40765 Egypt 1000 3000 61.3
## 15 Sovereignty uns… YES 102 Abyei 2011 3000 3.61
## # ℹ 3 more variables: Shape_Area <dbl>, Name_label <chr>,
## # geometry <MULTIPOLYGON [°]>Vamos a graficar el objeto countries:
La función names devolverá los nombres de todos los atributos almacenados.
Los objetos sf son objetos que cuentan con una columna especial de “geometría”:
st_geometry(countries)
## Geometry set for 15 features
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS: WGS 84
## First 5 geometries:
## MULTIPOLYGON (((30.47352 -1.057411, 30.47035 -1...
## MULTIPOLYGON (((39.37656 -4.721247, 39.37412 -4...
## MULTIPOLYGON (((33.02315 12.22673, 33.25586 12....
## MULTIPOLYGON (((40.21717 -10.29314, 40.21717 -1...
## MULTIPOLYGON (((34.09062 3.878785, 34.09114 3.8...Ya que las geometrías no contienen valorer únicos, se almacenan en una columna-lista, con cada elemento en la lista conteniendo la geometría de ese elemento. Las tres clases usadas para representar Simple Features son:
Si queremos graficar solo la geometría:
También podemos definir los límites para el eje x y el eje y del gráfico.
countries_geom <- st_geometry(countries)
plot(st_geometry(countries), axes = TRUE,
xlim = c(20, 40), ylim = c(1, 20))
Si queremos graficar un solo atributo:
Existen muchas más posibilidades: https://r-spatial.github.io/sf/articles/sf5.html
Cada fila contiene una única característica:
Para obtener una simple feature que contenga solo Etiopía:
ethiopia <- countries[7, ]
print(ethiopia)
## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 32.99994 ymin: 3.401109 xmax: 47.98618 ymax: 14.89416
## Geodetic CRS: WGS 84
## # A tibble: 1 × 10
## STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng Shape_Area
## <chr> <chr> <int> <chr> <int> <int> <dbl> <dbl>
## 1 Membe… NO 79 Ethiopia 1000 3000 50.4 92.9
## # ℹ 2 more variables: Name_label <chr>, geometry <MULTIPOLYGON [°]>
Para guardar la(s) característica(s) extraída(s) en una carpeta específica del sistema, podemos usar la función write_sf.
Puede guardar en otros formatos con st_write. Consulte el argumento driver en la ayuda.
Podemos deshacernos de la geometría con la función st_drop_geometry.
Calculando nuevos atributos (por ejemplo, mm/m2 al total por país):
Alternativamente, podemos usar dplyr:
countries %>%
mutate(precip_ml = runif(nrow(countries)),
precip_total = precip_ml * Shape_Area)
## Simple feature collection with 15 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 12.19545 ymin: -13.45568 xmax: 47.98618 ymax: 31.66734
## Geodetic CRS: WGS 84
## # A tibble: 15 × 13
## STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR Shape_Leng
## * <chr> <chr> <int> <chr> <int> <int> <dbl>
## 1 Member State NO 205 Rwanda 1000 3000 8.13
## 2 Member State NO 133 Kenya 1000 3000 48.5
## 3 Member State NO 74 South Su… 2011 3000 46.9
## 4 Member State NO 257 United R… 1000 3000 83.0
## 5 Member State NO 253 Uganda 1000 3000 23.2
## 6 Sovereignty uns… YES 61013 Ilemi tr… 1000 3000 2.84
## 7 Member State NO 79 Ethiopia 1000 3000 50.4
## 8 Member State NO 77 Eritrea 1000 3000 50.0
## 9 Member State NO 43 Burundi 1000 3000 9.12
## 10 Member State NO 68 Democrat… 1000 3000 99.0
## 11 Sovereignty uns… YES 40760 Hala'ib … 1000 3000 8.51
## 12 Sovereignty uns… YES 40762 Ma'tan a… 1000 3000 2.06
## 13 Member State NO 6 Sudan 2011 3000 81.9
## 14 Member State NO 40765 Egypt 1000 3000 61.3
## 15 Sovereignty uns… YES 102 Abyei 2011 3000 3.61
## # ℹ 6 more variables: Shape_Area <dbl>, Name_label <chr>,
## # geometry <MULTIPOLYGON [°]>, precip_mm <dbl>, precip_total <dbl>,
## # precip_ml <dbl>Cuando queremos manipular rásters y datos vectoriales simultaneamente, el paquete terra es muy útil. Para leer datos vectoriales con este paquete, usaremos la función vect.
Un ráster consiste en una matriz de celdas (o píxeles), donde cada celda contiene un valor o conjunto de valores que representan cierta información.
Los ráster con datos discretos a menudo utilizan un código para representar diferentes tipos de datos.
Para abrir un archivo ráster GeoTIFF (o similar), se utiliza el comando rast del paquete terra.
gleam.path <- "C:/User/L3_Rasters_and_Shapefiles_Data/GLEAM_annual/GLEAM_2009.tif"
gleam <- rast(gleam.path)Para poder realizar operaciones ráster, todos los archivos ráster involucrados deberán tener la misma geometría:
Los metadatos del archivo ráster se pueden visualizar escribiendo el nombre del objeto donde se almacenó (en este caso, gleam).
print(gleam)
## class : SpatRaster
## dimensions : 720, 1440, 1 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source : GLEAM_2009.tif
## name : GLEAM_2009
## min value : -23.93291
## max value : 2148.19751En R, la información sobre el ráster (sistema de coordenadas, extensión espacial, número de celdas, etc.) se puede obtener o introducir fácilmente usando la función crs del paquete raster:
crs(gleam)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"
crs(gleam) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"Para graficar un ráster usamos la función plot.
Alternativamente, podemos usar la función click para obtener el valor de un cierto píxel.
Para salir de la función click, presione Esc.
El paquete terra es capaz de trabajar eficientemente con un gran número de archivos ráster; por lo tanto, crear objetos de varias capas es extremadamente fácil.
El uso de objetos grillados de varias capas (stacks) es muy útil (por ejemplo, cuando se trabaja con rasters que involucran datos ambientales como precipitación, evapotranspiración y temperatura).
Queremos usar estos archivos anuales de ETa para:
Paso 1: Listar todos los archivos dentro de la carpeta.
ssebop.path <- "C:/User/Example/SSEBop_annual"
ssebop.files <- list.files(ssebop.path, full.names = TRUE, pattern = ".tif")head(ssebop.files)
## [1] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2003.tif"
## [2] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2004.tif"
## [3] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2005.tif"
## [4] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2006.tif"
## [5] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2007.tif"
## [6] "../Data/L3_Rasters_and_Shapefiles_Data/SSEBop_annual/SSEBop_2008.tif"Paso 2: Crear un stack a partir de todas las capas ráster individuales.
Paso 3: Calcular la media multianual.
Actualmente tenemos el ráster stack guardado como ssebop.
print(ssebop)
## class : SpatRaster
## dimensions : 8082, 7772, 16 (nrow, ncol, nlyr)
## resolution : 0.009652, 0.009652 (x, y)
## extent : -20.00845, 55.00689, -38.00535, 40.00211 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : SSEBop_2003.tif
## SSEBop_2004.tif
## SSEBop_2005.tif
## ... and 13 more source(s)
## names : SSEBop_2003, SSEBop_2004, SSEBop_2005, SSEBop_2006, SSEBop_2007, SSEBop_2008, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 3272, 3353, 3299, 3296, 3462, 3333, ...Al igual que con una lista, usamos dobles corchetes para acceder a un ráster específico (o varios rásters) desde un stack.
print(ssebop[[2]])
## class : SpatRaster
## dimensions : 8082, 7772, 1 (nrow, ncol, nlyr)
## resolution : 0.009652, 0.009652 (x, y)
## extent : -20.00845, 55.00689, -38.00535, 40.00211 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : SSEBop_2004.tif
## name : SSEBop_2004
## min value : 0
## max value : 3353Al igual que con una lista, usamos dobles corchetes para acceder a un ráster específico (o varios rásters) desde un stack.
print(ssebop[[1:3]])
## class : SpatRaster
## dimensions : 8082, 7772, 3 (nrow, ncol, nlyr)
## resolution : 0.009652, 0.009652 (x, y)
## extent : -20.00845, 55.00689, -38.00535, 40.00211 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : SSEBop_2003.tif
## SSEBop_2004.tif
## SSEBop_2005.tif
## names : SSEBop_2003, SSEBop_2004, SSEBop_2005
## min values : 0, 0, 0
## max values : 3272, 3353, 3299Recortar un ráster (función crop) extraerá la parte de un ráster dentro de un área rectangular.
Enmascarar un ráster (función mask), extraerá los píxeles dentro del shapefile.
La función crop extrae una porción de un conjunto de datos grillados según una entrada seleccionada (puede ser un archivo ráster o un archivo vectorial; es decir, SpatRaster o SpatVector, respectivamente).
Para recortar un ráster:
El parámetro snap se utiliza para alinear el archivo recortado con la extensión del elemento seleccionado que se utiliza para realizar el recorte.
Para enmascarar un ráster:
Para listar todos los valores únicos en un objeto ráster, se utiliza la función values.
-Como resultado se ontiene un vector
gleam_vals <- values(gleam.mask)
gleam_vals <- as.numeric(gleam_vals)
print(gleam_vals)
## [1] NA NA NA NA NA 823.3318 742.6441
## [8] 824.1913 NA NA NA 1063.5068 1031.8060 977.4069
## [15] 837.3028 793.3273 835.6875 NA NA 1141.6630 1031.7339
## [22] 938.8115 869.1260 827.8035 774.6710 820.8055 905.7816 NA
## [29] 1192.4106 1073.8735 914.8868 840.9797 838.4015 800.3232 828.6646
## [36] 908.3577 NA 1179.2241 1033.6238 878.3455 838.2838 795.3875
## [43] 769.3068 844.0502 860.8229 1094.2697 1071.3359 979.1405 876.9600
## [50] 837.6367 830.4340 897.6795 830.6163 840.5397 966.6479 974.6377
## [57] 1011.7329 880.3939 856.8268 NA NA NA NA
## [64] NA NA 957.3956 891.5947 873.2820 NA NA
## [71] NA NAPara realizar operaciones ráster, la geometría entre las capas debe ser la misma.
Se pueden utilizar diversos métodos.
Aquí nos enfocaremos en dos: interpolación bilineal e interpolación por el vecino más cercano (nearest neighbour).
La diapositiva anterior mostró el remuestreo a un origen diferente utilizando el mismo tamaño de celda de la grilla.
El remuestreo también se puede usar para cambiar la resolución espacial de la(s) capa(s).
Cómo remuestrear rásters en R. - Método nearest neighbour (rásters con variables categóricas o continuas)
print(land.cover)
## class : SpatRaster
## dimensions : 6002, 4802, 1 (nrow, ncol, nlyr)
## resolution : 0.008329863, 0.008330556 (x, y)
## extent : 10, 50, -15, 35 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : MCD12Q1_LC1_2009_001.tif
## name : MCD12Q1_LC1_2009_001
## min value : 1
## max value : 17Para remuestrear el ráster a la resolución espacial del objeto gleam, usaremos la función resample.
print(landcover.res)
## class : SpatRaster
## dimensions : 720, 1440, 1 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : MCD12Q1_LC1_2009_001
## min value : 2
## max value : 17
Cómo remuestrear rásters en R. Método bilineal (rásters con variables continuas).
print(chirps)
## class : SpatRaster
## dimensions : 2000, 7200, 1 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : chirps_monthly2009-01.tif
## name : chirps_monthly2009-01
## min value : 0.000
## max value : 1301.122
chirps.res <- resample(chirps, gleam, method = "bilinear")print(chirps.res)
## class : SpatRaster
## dimensions : 720, 1440, 1 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : chirps_monthly2009-01
## min value : 0.000
## max value : 1221.876
Después de haber modificado un rásters, es posible que deseemos exportarlo. Esto se puede hacer con la función writeRaster (parte del paquete terra).
x: objeto ráster. filename: nombre de archivo de salida (con la ruta de archivo completa). overwrite: argumento lógico.
Utilicemos la función writeRaster para guardar el ráster CHIRPS remuestreado del último ejercicio. Guardaremos el ráster como un archivo GeoTiff.
Pregunta: ¿Cómo podrías hacer que el nombre de la ruta de archivo cambie según el valor de una variable?
Una vez que los rásters tienen la misma geometría, es fácil realizar operaciones matemáticas.
Ahora, si tenemos estos rásters en un stack:
¿Cuál sería el resultado de las siguientes expresiones para la celda de cuadrícula resaltada:
Podemos identificar píxeles con ciertos atributos. Por ejemplo, las regiones con más de 200 mm de precipitación durante enero de 2009.
Para este propósito, estableceremos los valores mayores a 200 mm igual a 1, mientras que los valores menores o iguales a 200 mm se reemplazarán por un valor de 0.
Otro uso común de esta función es crear un ráster donde los valores son 1 o NA (de aquí en adelante, llamaremos a esto un ráster “unitario”). En este ejemplo, queremos calcular la precipitación promedio de las zónas áridas para enero de 2009 en los países del Nilo usando las zonas climáticas principales de Köppen-Geiger.
Ahora Encontraremos valores únicos usando la función values en combinación con la función unique.
| Valor | Clima |
|---|---|
| 1 | Tropical |
| 2 | Árido |
| 3 | Templado |
| 4 | Polar |
A continuación, establecemos los valores de la zona climática árida (valor = 2) a 1 y todo lo demás a NA.
print(arid)
## class : SpatRaster
## dimensions : 10548, 8367, 1 (nrow, ncol, nlyr)
## resolution : 0.004278282, 0.004278282 (x, y)
## extent : 12.1931, 47.98949, -13.45713, 31.67019 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : spat_vR6JxcbrnNc7Nbj_40545.tif
## name : NB_Major_CZ
## min value : 1
## max value : 1
print(chirps)
## class : SpatRaster
## dimensions : 2000, 7200, 1 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : chirps_monthly2009-01.tif
## name : chirps_monthly2009-01
## min value : 0.000
## max value : 1301.122Vamos a redimensionar el ráster unitario a la misma resolución espacial que el producto de precipitación.
arid <- resample(arid, chirps, method = "near")
print(arid)
## class : SpatRaster
## dimensions : 2000, 7200, 1 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : NB_Major_CZ
## min value : 1
## max value : 1Ahora, multiplicamos el ráster de precipitación por el ráster unitario de las zonas áridas y lo recortamos usando la capa climate.zones.
Finalmente, queremos obtener la precipitación media para la zona climática árida.
A veces, queremos extraer valores específicos (según las coordenadas) de un ráster. Para este propósito, utilizaremos la función extract del paquete raster. El primer paso es crear dos vectores con valores de latitud y longitud.
Ahora tenemos que convertir las coordenadas a un SpatVector.
Podemos graficar CHIRPS y las estaciones de la siguiente manera:
Para extraer los valores de las estaciones seleccionadas utilizaremos la función extract.
Como la función extract puede estar incluida en otro paquete, los dobles dos puntos (::) relacionan el paquete que se debe usar con la función respectiva. También podemos usar extract para extraer series temporales de un objeto stack.